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r~^ , The notion of duality between the hydrodynamic and kinetic (ghost) variables of lattice kinetic 

C^ • formulations of the Boltzmann equation is introduced. It is suggested that this notion can serve 

f^ ' as a guideline in the design of matrix versions of the lattice Boltzmann equation in a physically 

^N) , transparent and computationally efficient way. 
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I. INTRODUCTION 



In the last decade, the lattice Boltzmann (LB) method has developed into a very flexible and effective numerical 
technique for the simulation of a large variety of complex, fluid dynamical and non-equilibrium transport phenomena 
[l|. The LB method is based on a stream-and-coUide microscopic dynamics of fictitious particles, which stream with a 
T^ ' discrete set of velocities and interact according to local collision rules that drive the system towards a local equilibrium 
Sn . [ill 01 y| • Mathematically, this is formulated as the lattice Boltzmann equation (LBE) 

(1) 
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j/2 ' where /i(x. t) is the mean number of particles at position x and time i, moving along the lattice direction defined by 
the discrete velocity c^, (i, j — 1, ..., N). In the above. 









^ ' is a local equilibrium distribution, the discrete analogue of a Maxwellian distribution in continuum kinetic theory 
O truncated to second order in the mean flow velocity v, the Wi are a set of weights which satisfy ^^ wt — 1, and Cg is 
the speed of sound in the LB fluid. Greek indices denote Cartesian directions and the summation convention is implied. 
The low order velocity moments of the distribution function are related to the densities of mass, momentum and the 
deviatoric stress, {p, pva, Sap} = X^i /i{lj Cio,, Qia/?} where Sap + pc^^Sap = ^ap is the Eulerian momentum flux, and 
Qiap + c^s^ap = CiaCip. The higher moments of the distribution are related to the densities of rapidly relaxing kinetic 
degrees of freedom, variously called ghost or kinetic variables. Finally, Lij is a scattering matrix whose eigenvalues 
control the relaxation of the kinetic modes to their local equilibrium values. The null eigenvalues correspond to 
^jl . the eigenvectors associated with the conserved mass and momentum densities, while the leading non-zero eigenvalue 
associated with Qiap controls the viscosity of the LB fluid. 

Historically, the LBE in matrix form was derived as a Boltzmann approximation to the dynamics of lattice gas 
cellular automata 4]. It was then understood that the equilibria and collision matrix could be constructed indepen- 
f ^ \ dently of the underlying cellular automata microdynamics [5,] , and the lattice Boltzmann approach came into being. 
The collision matrix was reduced to the simplest possible form consistent with the macroscopic hydrodynamics in 
[a, 01, where the Bhatnagar-Gross-Krook (BGK) § form of the collision term was implemented on the lattice with 
Lij = T~^Sij . In the lattice BGK (LBGK) model, the fluid viscosity, which is the only transport parameter of interest, 
$H ' is given hy v — c^{t— 1/2). Even though it has always been clear that this simplification entails a crude approximation 
to the relaxation process (all modes relax at the same rate r), the LBGK equation, since its introduction, has held 
the mainstream in LB applications. In a parallel development, the collision matrix version of LBE [5, 9] has been 
revisited, optimized, and renamed the MTR (Multiple Time Relaxation) [rol, dl'l, in contrast to the single relaxation 
time implied by the LBGK equation. A number of authors have also made a strong case for the superiority of the 
MTR version over LBGK in terms of numerical accuracy and stability [13. [l3J. Yet, lattice BGK remains by far the 
most popular form of LBE to date. 

The limited of popularity of the MTR approach inspite of its superiority compared to the lattice BGK method, 
may be due to the lack of a general guiding criterion for the spectral decomposition of the collision matrix. In other 
words, it has not been clear a priori how to choose the eigenvectors of the matrix Lij which span the kinetic space 
of the discrete populations fi. This ambiguity arises because the conservation laws of mass and momentum fix only 
the hydrodynamic and transport subset of the eigenvectors, leaving the kinetic subset unspecified (see below). As 
a consequence, MTR models are dependent on both spatial dimension and the choice of the discrete velocity set c^. 



while the LBGK model is identical across both spatial dimension and choice of velocity set. Further, the notion of 
orthogonality of eigenvectors in the kinetic space can itself be defined in two distinct ways: the definition followed in 
[Sl using a weighted inner product, and that in [lO| using an unweighted inner product. For example, recent work 
on the shallow water equation 14], the fluctuating lattice Boltzmann equation [15], and on niultireflection boundary 
conditions use a set of eigenvectors which are orthogonal under the weighted inner product [16] . Clearly, it is important 
to understand how this non-uniqueness in the kinetic space arises and to provide a guiding principle in chosing the 
eigenvectors. In this work we shall propose such a guiding rule, by introducing the notion of duality in the kinetic 
space of the LB. 

In the next section we follow the notation in [l3| to highlight the structure of the kinetic space spanned by the 
eigenvectors and how a change of basis from populations to the moments reveals the dynamics of the various modes. 
We then introduce and illustrate the idea of duality with concrete examples. We show how a model in which all the 
ghost degrees of freedom are relaxed at the same rate [13, ll8[ may provide the best compromise between full MTR 
models where every mode has a separate relaxation time and the LBGK model. We end with a discussion on how the 
present work is relevant to the algorithmic improvement of the LB method. 

II. SPECTRAL REPRESENTATION OF THE COLLISION MATRIX 

A. Eigenvectors and eigenvalues 

For a general athermal DdQn LB model with n velocities in d space dimensions, the n x n collision matrix Lij 
has d + 1 null eigenvectors corresponding to the density and d components of the conserved momentum, d{d + l)/2 
eigenvectors corresponding to the stress modes, and n — (d + 1) — d{d + l)/2 eigenvectors corresponding to the 
ghost modes fl5J, [131 . The choice of the null and stress eigenvectors {l,Cia,Qiaf3} follows directly from the physical 
definition of the densities associated with them. Without specifying the exact analytical expression for the remaining 
eigenvectors, let us label a linearly independent set of the eigenvectors of the scattering matrix by {Af}, where 
a = 1 . . .n labels the eigenvector, and i = 1 . . .n labels the component of the eigenvector in along the i — th velocity 
direction. Then, we can define densities associated with the eigenvector Af as moments of the populations by 

V'"(x,t)=^/,(x,i)A^ (3) 

i 

For Af — {1, Cia^ Qiap} the densities are the mass, momentum and stress. The ghost eigenvectors are higher polyno- 
mials of the discrete velocities [9| . The discreteness of the kinetic space implies that, unlike in the continuum, only 
a finite number of polynomials can be linearly independent, being equal to the number of discrete velocities. For a 
model with n discrete velocities, the choice of the n linearly independent polynomials is thus not unique, but defined 
only upto a similarity transformation. Thus, the reason for the non- uniqueness in the spectral decomposition can be 
traced to the discreteness of the velocity space itself . Independent of the precise choice, the distribution function 
itself can be expanded in a linearly independent set eigenvectors which are polynomials of the discrete velocities 

M^,t)^w,J2^^iK,t)^ (4) 

a 

Consistency between the above two equations implies that the set of polynomials v4.^ are both orthogonal and complete, 

E, w^AfA'i^N'^S^", (5) 

Ea Af^VN-^S.,. (6) 

Crucially, with the definitions above [1^ , the eigenvectors Af form an orthogonal set under an inner product {A"' ,A'') = 
^iWiAfA'f. This inner product is identical to that introduced by Benzi et al Q, but distinct from the unweighted 
inner product {A°',A'') = "^^AfA^f used by d'Humieres and co-workers [10, [ill- The advantages of the present choice 
are discussed below. As indicated before, a useful categorisation of the polynomials consists of the d+1 polynomials 
{1, Cia} corresponding to the mass and momentum, the d{d -I- l)/2 quadratic polynomials Qia0 corresponding to the 
stress, and the remaining n — (d + I) — d(d + l)/2 cubic and higher order polynomials corresponding to the ghost 
variables. Correspondingly, the distribution function can be separated into contributions from the hydrodynamic, 
transport, and ghost moments 

/. = /f + /f + /f. (7) 



This motivates the introduction of projection operators [l7| which project the distribution function onto the hydro- 
dynamic, transport, and ghost subspaces, 

E p^fi = /f = «^^ E rAf/N- (10) 

j aGG 

(11) 

The exphcit form of the projection operators are 

i^f = w.,{l + c,-c,/cl) (12) 

Aj = WiQiapQjaff/'icj (13) 

P^'J = Y.wa^A^N'^ (14) 

aeG 

The discrete Maxwelhan is a nonhnear (quadratic) function of the distribution function, and thus Eq[T] is a only 
apparently linear, the nonlinearity being concealed in /?. An useful linearisation of the LB equation consists of 
neglecting the quadratic term in the discrete Maxwellian to yield a local equilibrium h^ which is linear in the mean 
velocity. 



'^? = -^(p+^)-E^f/. 



(15) 



In the linearised approximation for the equilibrium distribution, we have f^ ^ h^ = {P f)i and so the linearised 
LBE can now be written as, 

dth + c. • V/. = - E L,3 [fj - iP".f)j] = - E ^^f^ (16) 

j j 

where L^ = ^^ Lik{l — P^)kj is a right-projected collision matrix. In this form, it is clear that L^ by construction 
has eigenvectors of mass and momentum with zero eigenvalues. The form of the matrix, by itself, places no constraint 
on the eigenvalues of transport and ghost sectors. However, the requirmcnts of an extended range of hydrodynamic 
behaviour, stability and isotropy motivate an optimal construction of Lfj. As explained in the Introduction, the 
simplest possible model consists of a diagonal collision matrix L^ — Sij/r which implies that all the non-conserved 
modes relax at the same rate 1/r. This is the very popular LBGK approximation used in the literature. In the 
hydrodynamic regime, a scale separation exists between the relaxation of the conserved and non-conserved variables: 
the mass and momentum densities relax slowly, the stress and ghost variables relax rapidly. One variant of a model 
used by Ladd [ij] uses adjustable relaxation times for the stress modes, and identical unit relaxation times for the 
ghost modes, i.e. the ghosts are 'projected' out. One advantage of this approach is that the precise form of the ghost 
modes, which in general differ both in number and in form between LB models, need not be known. A generalisation 
of this model, with two relaxation times |J/7| reads, 

Lf^ = XP^^+aP,"^^ ail -P,f) + iX-a)P^^ (17) 

where the last follows from the completeness relations P^ + P^ + P'^ = 1. Since the precise form of the ghost 
projection operator, and hence the ghost eigenvectors is never needed in this formulation, it is clear that the linearised 
dynamics in this two-relaxation time model cannot depend on the precise choice of the ghost mode eigenvectors. 
The only way this model may be optimised is to tune the relaxation rate of the ghost modes in comparision to the 
stress modes. However, a model which allows separate relaxation times for each individual ghost mode has a greater 
flexibility and may be optimised to yield the best range of hydrodynamic behaviour |ll| . It needs careful analysis to 
see if the gain is enough to justify the loss of simplicity and generality that one obtains from the two relaxation model. 
Hydrodynamic behaviour is obtained when there are two propagating modes with a dispersion relation oj = Cgk+ivLk^ , 
ul = V + 3/2i^buifc being the longitudinal viscosity, and d — 1 diffusive modes with a dispersion relation bj = ivk^ . 
Both the speed of sound and the viscosities are assumed to be constant. 



B. Linear mode structure 

The hydrodynamic behaviour of the hnearised LBE is most conveniently analysed in the absence of boundaries when 
a Fourier mode decomposition is possible [i3i Hj [l3, [3 ■ It is important to note that the departure from hydrodynamic 
behaviour can arise from two distinct sources. The first is the choice of eigenvectors and relaxation times of the discrete 
velocity (but space and time continuous) LBE. This is the category of error arising from discretisation in velocity 
space. The second is that arising from the numerical integration of the LBE. This is the category of error arising 
from discretisation in space and time. The physical and numerical behaviour of the fully discretised LBE dynamics 
is a combination of both these sources of error. The present work, focussing as it does only on the kinetic space, has 
direct implications for errors arising out of discretisation of velocity space. The errors arising out of discretisation of 
space and time are relatively well understood from the numerical analysis of the hyperbolic differential equations. In 
particular, it is known that an Euler integration step of size At produces numerical diffusion, and thereby renormalises 
the viscosity to ly = c^,{t - At/2) 0|. 

To derive the dispersion relation we Fourier transform the linearised LBE to get 

9J. + zk.c./, = -^L,^/, (18) 

j 

At k = 0, the eigenmodes of the dynamics are the same as the eigenmodes of L^ . However, away from k = 0, neither 
the eigenmodes nor the eigenvalues are identical. For small k, an analytical expression for the eigenvalues may be 
obtained perturbatively [l7[ . For arbitrary k, a numerical solution is necessary. The dispersion relation is obtained 
by a Fourier transform in time, 

- iLo{k)n = - 5][*k • c,<5,, + L§f, (19) 

i 

Thus we need to obtain the eigenvectors and eigenvalues of the matrix 

My=zk-c,%+L^. (20) 

The dynamics in Eq llSI can equally well be written in terms of the densities using Eq|4] as 

b 

where matrix coupling the different modes is 

jsiaj^ab = ik • ^ w,A1A\c, (22) 

i 

It is worth noting that the hnearised LB dynamics can be written in either of the forms 

^t]^ = -^(^ + C),,/, (23) 

i 

dti^'' = -^(y^ + C)"V (24) 

b 

The dynamical equation in the fi basis diagonalises the advection operator Aij — ik • CiSij , while the dynamical 
equation in the V'° basis diagonalises the collison operator Cij = Lf-,. The eigenvectors of the dynamics are a 
combination of the fi and the ^". The dispersion relation equation can be conveniently non-dimensionalised by 
measuring time in units of the inverse of the relaxation rate for the stress modes r = A~^, and distance in units of 
CgT. The non-dimensionalised dispersion equation then takes the form 



zr!(q)/, = - ^[zq • ^% + Cf^f, (25) 



where fi = ujt is a non-dimensionalised frequency and q = kc^r is a non-dimensionalised wavevector. It should be 
noted that the non-dimensionalised collision matrix C^ = L^' /t now depends only on the ratio a/X of the relaxation 
rates of the ghost and stress eigenvectors. We shall use this non-dimensionalised form of the dispersion relation to 
obtain the numerical eigenspectrum of one of the LBE models presented below. 



III. DUALITY IN LATTICE KINETIC THEORY 

The symmetry principle of duality, which relates two different mathematical representations of the same physical 
theory, is a powerful tool in many areas of physical science. Duality is often use to map strongly interacting degrees 
of freedom to weakly interacting ones, thus facilitating an approximate, and often, even an exact solution of the 
problem. A celebrated example is the solution of Kramers and Wannier for the critical temperature of the Ising 
model [20|. To the best of our knowledge, dual symmetries do not appear to have played any major role in kinetic 
theory. In the context of the lattice Boltzmann schemes, we introduce duality not as an exact symmetry, but as 
a requirement on the structure of the kinetic space of the theory. Specifically, we require that the structure of the 
ghost subspace should mirror that of the hydrodynamic subspace, and consist of scalar densities and associated vector 
currents. In our notation, a LB kinetic space is exactly dual if each ghost field corresponds to a hydrodynamic field 
and a suitable transformation converts the ghost degrees of freedom into hydrodynamic degrees of freedom. If this 
exact correspondence is broken, but the ghost subspace still consists of sets of scalar densities and vector currents 
we say that the kinetic space is quasi-dual. The scalar densities and vector currents are taken to be even and odd 
functions of the discrete velocities respectively. Thus introduced, duality is a normative principle on the structure of 
the kinetic space of the LBE. The duality principle, as we show with several examples below, allows us to choose the 
eigenvectors of the collision matrix in a way which is both transparent and unique. 

A. Two dimensions 

Let us first consider the standard D2Q9 model with the usual set of velocities connecting the four nearest neighbours 
and the four next-nearest neighbours of the square lattice. Thus there are four velocities with unit modulus, another 
four with modulus two, which together with the zero velocity give the nine dynamical populations of the D2Q9 model. 
The kinetic space is spanned by eigenvectors corresponding to the mass and momentum, {A'^, Aj, Af} = {1^, Qj,, Ciy}. 
The next three natural eigenvectors associated with stress tensor are {A^,Af,A^} = {Qixx,Qixy,Qiyv} = {cL — 
c'^TCixCiyjcfy — Cg}. All of these are recognized as discrete velocity analogues of tensor Hermite polynomials |21|] . 
Without any physical considerations to guide us, the choice of three higher-order eigenvectors, associated with the 
ghost modes remains open. An obvious choice is the next series of tensor Hermite polynomials, that is QixxCix, 
QixxCiy, QiyyCix, QiyyCiy. It is immediately seen that due the identity cf^, = Cix, holding for the D2Q9 lattice, only 
two of these are linearly independent. This lack of linear independence, as we mentioned earlier, is due to the discrete 
nature of the velocities, giving identities like c|j, = Cix, which are absent in the continuum. To complete the kinetic 
space, one more eigenvector is required. It is immediately checked that, as a consequence of the D2Q9 identity 
cfa — c^ai o, — ^iUj out the fivc Hermite polynomials of order 4, only one is linearly independent, which we chose as 
QixxQiyy This then completes the construction of the remaining three ghost eigenvectors. 

In a very illuminating paper, Dellar [l4| proposes a different decomposition, based on the notion of ghost densities 
introduced in 9] . The first ghost eigenvector is of the form 

G° = 5. = (1, -2, -2, -2, -2, 4, 4, 4, 4) (26) 

and the remaining two are simply the corresponding 'currents', that is 

C^i = 9iCix, Gi = giCiy (27) 

The physical meaning of this choice is best highlighted by expressing gi in analytical form, that is 

„4 r 2 

where cf — cf^ -I- c^y. It is easily checked that the basis A° ... A^ ,A^ — G^ , A'' — G^,A^ = G'^ is orthogonal under the 
weighted scalar product [A"- , A^) = ^^ WiA°:A\, where Wq = 4/9, Wi_4 = 1/9 and Ws-g = 1/36 are the usual D2Q9 
weights. It is also to be noted that, owing to the D2Q9 identities, the ghost eigenbasis can also be written as 

(29) 
(30) 
(31) 
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Surprisingly, then, G^ = gi is a. fourth order lattice Hermite polynomial, while Gj = giCix and Gf 
of being fifth order lattice Hermite polynomials turn out to be third order lattice Hermite polynomials. This fulfils 
exactly the duality principle introduce above: the kinetic space is decomposed into a set of eigenvectors corresponding 
to conserved, transport and ghost moments; the ghost degrees of freedom correspond to an even scalar density and 
two odd vector currents and are in one-to-one correspondence with the hydrodynamic degrees of freedom; and as we 
show below, the ghost and hydrodynamic degrees of freedom are related by a suitable transformation. 

The duality in the decomposition is beautifully illustrated by the diamond structure of the D2Q9 eigenvectors 
shown in Table 1. The density and the two momenta are matched by a ghost density and two ghost currents. The 
dynamical behaviour of these degrees of freedom are of course quite different, as is revealed by displaying the LBE 
dynamics in the basis of moments. The kinetic moments associated with the present choice of eigenvectors is 

{p,pVa,Sa[3,p',j'a} == ^ /i{l, Qa , Qia/S, 5i, 5iCia } (32) 

i 

The primed quantities correspond to ghost density and its currents. The decompostion of fi as the sum of a hydro- 
dynamic, transport and ghost components is 

/f = w,{p + t^) (33) 

fi = M TTT-^) (34) 



2cf 
■^w^g^ip + —^ 



f? = 7^^9^ip' + '-^) (35) 



From the above expressions it is clear that, to within a scale factor, the ghost sector is transformed into the conserved 
sector under the duality transformation 1^ ^^ gi- 

To physically interpret the above decomposition, we first note that the combination w'^ = Wigi may be interpreted 
as the weight associated with the ghost degrees of freedom. Then, the weights of the hydrodynamic modes sum to 
unity "^Wi = 1, while the weights of the ghost modes sum to zero X^i ''^i ~ 0. This last result combined with the 
fact that ghost density is even in the velocities gi = gi* , where Ci* ~ — Ci, indicates that the ghosts correspond to 
oscillatory eigenvectors familiar in quantum and statistical mechanics, where they represent excitations above the 
ground state or above equilibrium. The ghost degrees of freedom are thus non-equilibrium excitations carried by 
even, oscillatory eigenvectors. The even and odd character of eigenvectors can be exploited to classify the entire set 
of moments into two categories: even moments representing densities, and odd moments representing currents. Odd 
moments, representing currents, vanish at global equilibrium by symmetry. The even moments are not constrained 
to vanish by symmetry arguments. However, since the kinetic modes have no projection onto the global equilibrium 
distribution function f^ = Wip, they can be conveniently chosen to vanish at equilibrium. This is one of the principal 
advantages of using a set of eigenvectors which are orthogonal under the weighted inner product. 

By interpreting Wi as 'masses' of the hydrodynamic modes, the wl can be identified with 'masses' of ghost modes. 
By construction, since they sum up to zero, some of these masses ought to be negative. For instance, the ghost density 
can be rewritten as an alternating sum of the populations associated with the three energy levels c^ = 0, 1, 2, that is 
p' ^ fo- 2(/i + /2 + /a + fi) + 4(/5 + /e + /? + /s) = Po - 2pi + 4/92, where j = 0, 1, 2 refer to the j-th energy level. 
Being the sum of populations, each of the partial densities is strictly non-negative at all times, but the combination 
of alternating coefficients is a signed quantity, p' , which is zero only at equilibrium. The duality is made even more 
apparent by defining the reduced distribution function (pi = fi/wi, thus writing p — ^^ Wi<pi and p' = ^^ 'w'^(|)i. Since 
Wi is the lattice analogue of global equilibrium distribution, w[ may also be interpreted as a measure of the global 
departure from equilibrium. 

The ghost currents are a measure of the skewness of the kinetic distribution function, which is non zero only out 
of equilibrium. Being based on this equilibrium ^^ non-equilibrium duality, the ghost decomposition shows that the 
higher-order excitations, keeping the system away from equilibrium, can be structured exactly like their hydrodynamic 
counterparts. It should be appreciated the duality is structural and not dynamical: it is broken at various levels, 
starting with the prefactors defining the ghost density and current, because the norm of the dual hydrodynamic versus 
ghost eigenvectors is not the same. In particular, this implies that the ghost kinetic tensor is not isotropic, as one 
can easily check by a direct calculation: P^^ = AP^x and P' = ^AP^y This is not surprising, since equilibrium and 
non-equilibrium are not physically equivalent. 

However, a dynamical transformation in time, A ^^ 1/A turns perfectly conserved modes (infinite-lifetime) into 
perfectly non-conserved ones (zero-lifetime). What this means is that the distinction between equilibrium and non- 
equilibrium modes is not dictated by the structure of the kinetic space, but only by the actual values of the lifetimes 
of the excitations supported by this equation. In this respect, we expect a signature of this dynamical duality in the 
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TABLE I: The arrangement of the eigenvectors of the D2Q9 lattice in a diamond structure. There is an exact duahty about 
the transport sector (T) with the conserved hydrodynamic (C) and ghost (G) sectors, transforming into each other under the 
interchange of weights (see text). 
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TABLE II: The hierarchical tree of moments for the 16 speed dual lattice D{3) defined in the text. The first four levels are 
constructed using the usual lattice Hermite polynomials. The remaining three levels are completed using the duality prescription 
starting with an even, sixth-order scalar density 

form of a mirror symmetry A/cr <-> a/X in the dispersion relation for the two-relaxation time model introduced earlier. 
Numerical dispersion relations presented in the next section do show evidence of such a symmetry. 

The structural duality of the kinetic space is broken dynamically by the different eigenvalues assigned to the 
hydrodynamic and ghost modes. The hydrodynamic sector sustains itself even in the absence of ghost (standard 
macroscopic hydrodynamics), whereas the ghosts, because of finite lifetimes assigned to them, do not survive without 
a forcing from the hydrodynamic modes. Indeed, in the absence of the hydrodynamic feedback on P^a, the ghost sector 
would rapidly estinguish, because neither its density nor its current are conserved in time. Of course, such dynamical 
asymmetry can always be removed by choosing the ghost eigenvalues equal to be zero. In fact, there is even some 
evidence that long-lived ghosts may prove beneficial to the numerical stability of short-scale hydrodynamics in fluid 
turbulence [22] . This is a prescription to keep the Boltzmann distribution away from equilibrium for an indefinitely 
long time, leading to anomalous relaxation. This prescription clearly violates the normal ordering between slow, 
hydrodymic and fast, kinetic modes, as it corresponds to enforcing additional conservation laws with no counterpart 
in the real molecular world. Hence, such a procedure can only be justified as an effective interaction between collective 
degrees of freedom, as for example in lattice kinetic equations for turbulent flows. 

B. Higher order lattices in two dimensions 

In two dimensions, lattices with more than nine velocities are used in thermal LBE and in applications in microflu- 
idics and multiphase flow. Duality can be used to generate an optimal kinetic space for these higher order lattices 
as well. Let us denote by D{s) the lattice corresponding to a hierarchical tree of eigenvectors, with 2s -t- 1 levels and 
symmetric about the s + 1th level. (see Table 2 and 3). Clearly, this hierarchy contains (s + 1)^ independent moments. 
With this deflnition, D{Q) corresponds to the lattice with a single zero-speed rest particle, -0(1) to the D2Q4: lattice, 
and D{2) with the D2Q9 lattice. -0(3), which represents a higher order lattice in the present terminology, consists 
of 16 velocities with four velocities each of modulus 1,2,4 and 8. For this rather complicated lattice, the duality 
prescription proceeds by flrst chosing the usual eigenvectors corresponding to the conserved and transport sectors. 
Proceeding to the 3rd level, the eigenvectors of the type ciaQip-y (and permutations) now turn out to be linearly inde- 
pendent. The remaining moments are constructed in a top-down fashion, beginning with a sixth-order scalar density 
PQi = Ac^ — Bcf + Ccj — li from which two currents peiCia and three tensorial densities psiQixx, PeiQixy, PeiQiyy 
can be used to complete the hierarchy. The expansion coefficients A, B, C for the scalar density can be computed by 
requiring orthogonality to the lower eigenvectors. 

The next member of the hierarchy is -D(4), corresponding to the lattice of 25 speeds which has recently been 
shown to have 8th order isotropy in its spatial behaviour. As with D{3), the eigenvectors upto the 3rd level are the 
conserved, transport, and tensorial Hermite polynomials CiaQip-y and permutations. Again, on the 25 velocity lattice, 
the permutations give rise to independent eigenvectors. The remaining eigenvectors constructed bottom-up starting 
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TABLE III; The hierarchical tree of moments for the 25 speed dual lattice -D(4) defined in the text. The first four levels are 
constructed using the usual lattice Hermite polynomials. The remaining three levels are completed using the duality prescription 
starting with an even, eighth-order scalar density 

from an even, scalar eighth-order density psi = Ac^ — Bc^ + Ccf — Dcf + 1^, from which two currents psiCix, s psiCiy, 
three tensorial densities psiQixx, PsiQixy, PSiQiy, and four tensorial currents pstUxxxx, PsiUxxyy, PsiUxyyy, PsiUyyyy 
can derived, thus completing the list of 25 independent eigenvectors. 

Both of the above examples show that the duality prescription offers a transparent method of choosing and ordering 
the set of eigenvectors in lattices which are more complicated than the most commonly used D2Q9 lattice in two 
dimensions. It has recently been shown that the 16 and 25 speed lattice with a proper choice of weights, provide 6th 
and 8th order isotropy respectively. [23, [2J| . Since the choice of weights is intimately related to both weighted inner 
product and the duality prescription, it is possible that there is fundamental link between isotropy and the duality 
prescription. 

C. Three dimensions 

In three dimensions, the model which is potentially exactly dual is the D3Q1A model, which has the usual 4 
hydrodynamic degrees of freedom, 6 transport degrees of freedom, leaving 4 ghost degrees of freedom to match the 
hydrodynamic ones. However, since the Z53Q14 model is not used in practice, we pass on instead to the analysis of the 
most common D3Q19 model. Here, of course, the kinetic space can only be quasi-dual, since there are 9 ghost degrees 
of freedom. To choose them according to the duality prescription, several possibilities can be explored. With one ghost 
density quartic in the velocity, and three associated currents which are quintic, we obtain 4 ghost eigenvectors, leaving 
5 free. This allows three independent components of the ghost momentum-flux tensor {xy,xz,yz), plus another two, 
which must necessarily come from a higher Hermite level. This seems to be a rather obscure and unpromising avenue. 
A better possibility is to select two ghost densities along with their currents, leaving the third one 'naked', i.e. without 
independent degrees of freedom for the current. Retaining three quartic ghost densities with their respective currents 
is unviable, for it gives a total of twelve eigenvectors, three too many. From these considerations, it appears that 
the 2 (dressed) plus 1 (naked) density representation comes closest to fulfillingl the duality program. It is interesting 
to note that, apart from the third, naked, density, this is precisely the early decomposition adopted in [9'], based on 
the 3d projection of 4d face-centered hypercube (24 speeds in d = 4, 18 in d = 3). The explicit form of the chosen 
eigenvectors is given in the Appendix. 

IV. NUMERICAL RESULTS 

We now present a numerical calculation of the dispersion relation of the two-relaxation time lattice Boltzmann model 
with the duality-prescribed choice of eigenvectors. The dispersion relation is obtained by numerically computing the 
eigenvalues and eigenvectors of the matrix My. As explained previously, with a suitable rescaling, the only parameter 
in the problem is a/A, the ratio of the relaxation rates of the kinetic and stress degrees of freedom. For a = X 
the collision term reduces to the LBGK diagonal collision operator. The imaginary parts of the eigenvalues for the 
D2Q9 model are shown in Fig.l , clearly showing the presence of hydrodynamic modes (relaxation rates vanish as 
wavenumber goes to zero) as well as non-hydrodynamic modes (relaxation rates remain finite as wavenumber goes to 
zero) . The scale separation between the relaxation rates of the hydrodynamic and non-hydrodynamic modes becomes 
progressively smaller with increasing wavenumber, and there is considerable overlap at around g = tt. This is fairly 
plausible, since g = tt is the value at which the wavelength becomes comparable with the mean free path CgT, so that 
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FIG. 1: Dispersion relation of the D2Q9 matrix LBE for r = a — 1. This is identical to the BGK model. Note the large 
overlap of the hydrodynamic (lower curves) and non-hydrodynamic (upper curves) modes. 




FIG. 2: Dispersion relation of the D2Q9 matrix LBE for r = 
the stress modes, leading to poor hydrodynamic behaviour. 



1, a — 2. The ghost modes are forcibly made to relax slower than 



the distinction between hydrodynamics and kinetic modes fades away. This lack of scale separation is responsible for 
the poor range of hydrodynamic behaviour of the LBGK models, a fact that was correctly noted earlier ^^. With 
a = ^A, the overlap between the hydrodynamic and non-hydrodynamic modes in Fig. 2 is even greater, indicating a 
further reduced range of hydrodynamic behavior, compared to the LBGK models. On the other hand, for a = 2A, 
we see in Fig. 3 a clean separation between the hydrodynamic and kinetic degrees of freedom, and it is in this range 
of parameters that we expect the best hydrodynamic behavior of the two-relaxation time LB model. An interesting 
qualitative feature that emerges from comparing Fig. 2 and Fig. 3, is that the eigcn-frequencies are almost 'dual' to 
each other, in the sense that the dispersion curves are approximately identical after a reflection about the ordinate 
and a rescaling by — . The structural duality reflects itself in the dynamical behaviour if the relaxation times are 
chosen appropriately. 

This supports our earlier assertion that the structure of the hydrodynamic and non-hydrodynamic modes are dual 
to each other. It also illuminates the physical behavior of ghost modes, whose dynamics appears to be characterized 
by a competition between global decay as characterised by a, and instabilities driven by negative diffusion as indicated 
by the negative curvature of the ghost dispersion relations. Thus the ghost modes decay globally, but driven by the 
negative diffusion, concentrate around thinner and thinner regions of space, and can thereby undermine the high 
frequency high wavenumber stability of the system. Good hydrodynamic behavior is thus expected whenever global 
decay proceeds sufficiently fast to deplete the ghost energy before this energy has time to cascade to high frequencies. 
This picture suggests a number of interesting questions for future studies. First, it would be interesting to explore 
whether the use entropic methods [2^ simply accelerates the ghost decay, or rather turns ghosts into stable modes. 
Second, following [221, i^ would be interesting to study whether the dual-decomposition can help designing the ghost 
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FIG. 3: Dispersion relation of the D2Q9 matrix LBE for t — 1, a = 1/2. The ghost modes relax twice as fast the stress modes. 
There is a clean seperation of time scales and enhanced hydrodynamic behaviour compared to the BGK model. 

dynamics in such a way as to absorb energy bursts from the hydrodynamic component, as they occur in a turbulent flow 
(intermittency). This could be achieved, for instance, by promoting ghost eigenvalues to dynamical fields responding 
self-consistently to the local dynamics of the turbulent flow, as it is currently done with the transport eigenvalues r 
in the kinetic modeling of fluid turbulence [26| . Finally, we note that the long- wavelength dynamics obtained with 
the present choice of eigenvectors is, by construction, isotropic at order k'^ and Galilean invariant. 

V. CONCLUSION 

In this paper we have developed the notion of duality between the hydrodynamic and ghost sectors of lattice kinetic 
equations, as a guiding criteria to resolve the ambiguities which arise in the practical construction of LB models in 
matrix form. Our main prescription is that the ghost sector should be constructed, in analogy with the hydrodynamic 
sector, to consist of density-current pairs. This prescription is exactly realised in the D2Q9 model, where in addition, 
the ghost and hydrodynamic sectors can be interchanged by a suitable swapping of weights. For higher order lattice in 
and in higher dimensions the kinetic degrees of freedom are more numerous then the hydrodynamic ones thereby ruling 
out an exact correspondence between the two. However, the duality prescription still provides an useful ordering of 
the eigenvectors into a quasi-dual kinetic space. The duality principle presented in this paper has been used previously 
in constructing the kinetic space of the fluctuating lattice Boltzmann equation [l5|. It has also been recently used 
to compare the accuracy of multireflection boundary conditions with both weighted and un-weighted eigenvectors 
[16| . We hope the duality principle as introduced here will provide an impetus to further developments in the matrix 
formulation of the lattice Boltzmann method. 



VI. APPENDIX 



For easy reference we present the eigenvectors of the D2Q9 and D3Q19 models chosen according to the duality 
prescription with weighted inner product. The table is arranged according to the conserved (C), transport (T) and 
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